#include <bits/stl_alloc.h>
#include <stdexcept>
#include <string>
#include <algorithm>
#include <iostream>
#include <utility>
#include <vector>
#include <sstream>
#include <fstream>
#include <deque>
#include <boost/lexical_cast.hpp>
#include <boost/tuple/tuple.hpp>

# define DEBUG_LEVEL 9
#include <debug_func.h>
#include <perlutil.h>
#include <bl2seq_mine.h>
#include <portable_timer.h>
#include <batch_job.h>
#include <progress_indicator.h>
#include <pg_copy_guard.h>
#include <pg_db_cursor.h>
#include <bcpp_parsefasta.h>
#include <bcpp_read_seq.h>
#include <stlhelper.h>
#include <string_hash.h>
#include <print_error.h>
#include <ext/hash_map>
#include <gc.h>

#include "get_arg.h"

using BioCPP::t_collect_sequences;
using BioCPP::t_collect_sequences_report_errors_to_strm;
using BioCPP::eSEQ_NODUPLICATES;
using BioCPP::eSEQ_DUPLICATES;
using std::cerr;
using std::cout;
using std::string;
using std::vector;
using std::deque;
using std::istream;
using std::ifstream;
using std::sort;
using std::make_pair;
using std::pair;
using std::ostringstream;
using boost::lexical_cast;
using boost::tie;




//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		Main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


	
	

//________________________________________________________________________________________

//	t_collect_sequences_from_strm

//________________________________________________________________________________________
struct t_dna_freq_from_strm
{
	t_progress_indicator&				dots;
	BioCPP::t_valid_char				valid_seq_char;
	t_dna_freq_from_strm(	t_progress_indicator&	dots_,
							const std::string&		valid_sequence_chars)
			:
			dots(dots_),
			valid_seq_char(valid_sequence_chars),
		 {	}



	// helper function
	static void filter_sequence(std::string& sequence, const BioCPP::t_valid_char& valid)
	{
		sequence.erase(transform_if(sequence.begin(), sequence.end(), sequence.begin(),
				BioCPP::t_sequence_tolower(), valid), sequence.end());
	}
	
	// callback
	bool operator()(std::string& accession, std::string& sequence,
						unsigned acc_filepos,
						unsigned seq_filepos);
};

bool t_dna_freq_from_strm::operator()(std::string& accession, std::string& sequence,
											unsigned acc_filepos, unsigned seq_filepos)
{
	// print progress dots
	++dots;

	// get line counts for error messages
	unsigned curr_acc_beg_line = acc_beg_line;
	acc_beg_line += std::count(sequence.begin(), sequence.end(), '\n') + 1;

	// remove invalid characters
	filter_sequence(sequence, valid_seq_char);

	// bad sequence!!!
	if (!sequence.length())
	{
		std_print_error("Empty sequence at line #" + ultoa(curr_acc_beg_line));
		return true;
	}
	if (!accession.length())
	{
		std_print_error("Empty accession at line #" + ultoa(curr_acc_beg_line));
		return true;
	}



	return true;
}

int main(int argc, char* argv[])
{
	debug_func_cerr( 8 );
	
	t_program_args args;

	try
	{

	// don't use C streams
	std::ios::sync_with_stdio(false);

	// get command line arguments
	if (!process_command_line_options(argc, argv, args))
		return 1;

	

	//
	// read sequences
	//
	t_collect_sequences seq_db(args.align_mode == eALIGN_PAIRWISE ? eSEQ_DUPLICATES : eSEQ_NODUPLICATES);
	{
		args.VERBOSE() << "\tReading sequences from STDIN...\n";
		t_progress_indicator dots(args.VERBOSE(), 10);
		t_collect_sequences_report_errors_to_strm report_errors(args.CERR());
		seq_db.init_from_stream(std::cin, args.regex, dots, report_errors);
		// put into blast format
		seq_db.bl_prepare_all_sequences();
	}
	args.VERBOSE() << "\tDone...\n";
	args.VERBOSE() << "\t" << seq_db.size() << " sequences\n";
	
	
	vector<pair<unsigned, unsigned> > 	sequence_pairs;
	//
	// read sequences pairs
	//
	if (args.align_mode == eALIGN_SPECIFIED)
	{
		assert(args.file_pairs.length());
		args.VERBOSE() << "\tReading sequence pairs...\n";
		ifstream f_seq_names(args.file_pairs.c_str());
		if (!f_seq_names)
		{
			args.CERR() << "Error:\n\tCould not open " + args.file_pairs + " to retrieve "
							" file pairs\n";
			return 1;
		}
		read_seqid_pairs(f_seq_names, seq_db, sequence_pairs, args.CERR());
	}
	
	//
	// all vs all
	//
	else if (args.align_mode == eALIGN_ALLVSALL)
	{
		args.VERBOSE() << "\tAll vs all comparisons...\n";
		unsigned sz = seq_db.size();
		for (unsigned i = 0; i < sz - 1; ++i)
			for (unsigned j = i + 1; j < sz; ++j)
				sequence_pairs.push_back(make_pair(i, j));
	}

	//
	// sequence pairs not specified: successive pairs
	//
	else
		//args.align_mode == eALIGN_PAIRWISE
	{
		args.VERBOSE() << "\tcomparing successive pairs...\n";
		unsigned sz = seq_db.size();
		for (unsigned i = 0; i < sz - 1; i += 2)
			sequence_pairs.push_back(make_pair(i, i+1));
	}

	//
	//	blast options
	//
	t_portable_timer timer;
	t_blast_options blast_options(args.seg, args.cutoff, args.blosum);
	
	
	args.VERBOSE() << "\tAligning " << sequence_pairs.size() << " sequence pairs ...\n";
	
	
	if (args.results_header)
		cout	<< "seq1\tseq2\t"
				<< t_blast_results::get_tab_delimited_column_names() << "\n"
				<< "----\t----\t"
				<< t_blast_results::get_tab_delimited_column_name_underscores() << "\n";

	unsigned cnt_results = 0;
	{
		t_progress_indicator dots(args.VERBOSE(), std::min(20U, sequence_pairs.size() / 150),
																	sequence_pairs.size());
		dots.use_timer();
	
		//__________________8______________________________________________________________
		// 		align all genes in pairs
		//________________________________________________________________________________
		for (unsigned i = 0; i < sequence_pairs.size(); ++i)
		{
			unsigned seq1 = sequence_pairs[i].first;
			unsigned seq2 = sequence_pairs[i].second;
			if (align_sequences(dots, blast_options, args.lencutoff, seq_db, seq1, seq2, args))
				cnt_results++;
		}
	}
	
	args.VERBOSE() << cnt_results << " sequence pairs out of "
					<< sequence_pairs.size() << " could be aligned.\n";
	


	}
	catch (std::exception& e)
	{
		std_print_error(args.CERR(), e.what());
		return 1;
	}
	catch(...)
	{
		return 1;
	}

	return 0;
}

